
import os
import sys
import numpy as np
import pandas as pd
import sklearn
from sklearn.tree import DecisionTreeClassifier
from sklearn.ensemble import RandomForestClassifier, RandomForestRegressor, AdaBoostClassifier
from sklearn.svm import SVC, SVR
from sklearn.model_selection import train_test_split, StratifiedShuffleSplit, cross_val_score, GridSearchCV
from sklearn.metrics import confusion_matrix, ConfusionMatrixDisplay, roc_auc_score, roc_curve, RocCurveDisplay, mean_squared_error, mean_absolute_error, auc, f1_score
from sklearn.preprocessing import StandardScaler, label_binarize, LabelBinarizer
from sklearn.linear_model import LogisticRegression
from sklearn.linear_model import LinearRegression, Ridge, Lasso, RidgeClassifier
import scipy
import scipy.stats as stats
import statsmodels.api as sm
import matplotlib
import matplotlib.pyplot as plt
from matplotlib.colors import Normalize
import seaborn as sns


print(f"Python version: {sys.version}")

packages = [np, pd, sklearn, scipy, sm, matplotlib, sns]
for package in packages:
    print(f"{str(package).partition('from')[0]} using version: {package.__version__}")

Python version: 3.11.7 (main, Dec  4 2023, 18:10:11) [Clang 15.0.0 (clang-1500.1.0.2.5)]
<module 'numpy'  using version: 1.26.3
<module 'pandas'  using version: 2.1.4
<module 'sklearn'  using version: 1.3.2
<module 'scipy'  using version: 1.11.4
<module 'statsmodels.api'  using version: 0.14.1
<module 'matplotlib'  using version: 3.8.2
<module 'seaborn'  using version: 0.13.2


# Import data
current_wdir = os.getcwd()
data_folder = current_wdir + '/Data'
sleep_data_folder = current_wdir + '/Data/Sleep_Data'

# Iterate through each file in .Data/ and add it to a dataframe.
file_path = [f'{data_folder}/{file}' for file in os.listdir(data_folder) if '.csv' in file]
# Using a list to concat the dfs with index_col allows to easily merge based on 'day'. More memory usage but fine for this project.
biometric_df = pd.concat([pd.read_csv(file, index_col = 'day') for file in file_path], join = 'outer', ignore_index = False, axis = 1).reset_index()

# Iterate through each file in .Data/Sleep_Data and add it to a dataframe.
# Separated sleep data as it potentially has multiple entries per day. Will merge them later.
file_path_sleep = [f'{sleep_data_folder}/{file}' for file in os.listdir(sleep_data_folder) if '.csv' in file and 'daily' not in file]
sleep_df = pd.concat(map(pd.read_csv, file_path_sleep), join = 'outer', ignore_index = False, axis = 1)

# Import ground truth label sleep score.
file_path_daily_sleep_score = [f'{sleep_data_folder}/{file}' for file in os.listdir(sleep_data_folder) if '.csv' in file and 'daily' in file]
daily_sleep_score_df = pd.read_csv(file_path_daily_sleep_score[0])
daily_sleep_score_df = daily_sleep_score_df[['score', 'day']]


biometric_df.dtypes.value_counts()

int64      25
float64    13
object      4
Name: count, dtype: int64


sleep_df.dtypes.value_counts()

float64    33
object      9
int64       5
Name: count, dtype: int64


# First detailed looks
print('biometric_df info:\n------------------\n', biometric_df.shape)
display(biometric_df.head(2))
print(biometric_df.columns)
print('\n\nsleep_df info:\n------------------\n', sleep_df.shape)
display(sleep_df.head(4))
print(sleep_df.columns)

biometric_df info:
------------------
 (369, 42)

Index(['day', 'active_calories', 'average_met_minutes',
       'equivalent_walking_distance', 'high_activity_met_minutes',
       'high_activity_time', 'inactivity_alerts', 'low_activity_met_minutes',
       'low_activity_time', 'medium_activity_met_minutes',
       'medium_activity_time', 'meters_to_target', 'non_wear_time',
       'resting_time', 'sedentary_met_minutes', 'sedentary_time', 'steps',
       'target_calories', 'target_meters', 'total_calories', 'score',
       'class_5_min', 'contributors_meet_daily_targets',
       'contributors_move_every_hour', 'contributors_recovery_time',
       'contributors_stay_active', 'contributors_training_frequency',
       'contributors_training_volume', 'met_1_min', 'ring_met_1_min', 'score',
       'temperature_deviation', 'temperature_trend_deviation',
       'contributors_activity_balance', 'contributors_hrv_balance',
       'contributors_previous_day_activity', 'contributors_previous_night',
       'contributors_recovery_index', 'contributors_resting_heart_rate',
       'contributors_sleep_balance', 'contributors_body_temperature',
       'spo2_percentage'],
      dtype='object')


sleep_df info:
------------------
 (730, 47)

Index(['average_breath', 'average_heart_rate', 'average_hrv', 'awake_time',
       'bedtime_end', 'bedtime_start', 'day', 'deep_sleep_duration',
       'efficiency', 'latency', 'light_sleep_duration', 'lowest_heart_rate',
       'movement_30_sec', 'period', 'rem_sleep_duration', 'restless_periods',
       'score', 'segment_state', 'sleep_midpoint', 'time_in_bed',
       'total_sleep_duration', 'type', 'sleep_phase_5_min', 'timezone',
       'bedtime_start_delta', 'bedtime_end_delta', 'midpoint_at_delta',
       'heart_rate_5_min', 'hrv_5_min', 'contributors_total_sleep',
       'contributors_deep_sleep', 'contributors_rem_sleep',
       'contributors_efficiency', 'contributors_latency',
       'contributors_restfulness', 'contributors_timing',
       'readiness_contributors_activity_balance',
       'readiness_contributors_body_temperature',
       'readiness_contributors_hrv_balance',
       'readiness_contributors_previous_day_activity',
       'readiness_contributors_previous_night',
       'readiness_contributors_recovery_index',
       'readiness_contributors_resting_heart_rate',
       'readiness_contributors_sleep_balance', 'readiness_score',
       'readiness_temperature_deviation',
       'readiness_temperature_trend_deviation'],
      dtype='object')


print('\n\ndaily_sleep_score_df info:\n--------------------------\n',daily_sleep_score_df.shape)
display(daily_sleep_score_df.head(4))


daily_sleep_score_df info:
--------------------------
 (355, 2)


# Confirm there are no duplicates of days in the biometric or sleep score data.
print(biometric_df['day'].duplicated().sum())
print(daily_sleep_score_df['day'].duplicated().sum())

0
0


# Remove contributor and score columns (multiple scores involved, will add truth label later).
biometric_df = biometric_df.loc[:, ~biometric_df.columns.str.contains('contributors|score')]
sleep_df = sleep_df.loc[:, ~sleep_df.columns.str.contains('contributors|score')]

print(biometric_df.shape)
print(sleep_df.shape)

(369, 26)
(730, 30)


print(biometric_df.columns, sleep_df.columns)

Index(['day', 'active_calories', 'average_met_minutes',
       'equivalent_walking_distance', 'high_activity_met_minutes',
       'high_activity_time', 'inactivity_alerts', 'low_activity_met_minutes',
       'low_activity_time', 'medium_activity_met_minutes',
       'medium_activity_time', 'meters_to_target', 'non_wear_time',
       'resting_time', 'sedentary_met_minutes', 'sedentary_time', 'steps',
       'target_calories', 'target_meters', 'total_calories', 'class_5_min',
       'met_1_min', 'ring_met_1_min', 'temperature_deviation',
       'temperature_trend_deviation', 'spo2_percentage'],
      dtype='object') Index(['average_breath', 'average_heart_rate', 'average_hrv', 'awake_time',
       'bedtime_end', 'bedtime_start', 'day', 'deep_sleep_duration',
       'efficiency', 'latency', 'light_sleep_duration', 'lowest_heart_rate',
       'movement_30_sec', 'period', 'rem_sleep_duration', 'restless_periods',
       'segment_state', 'sleep_midpoint', 'time_in_bed',
       'total_sleep_duration', 'type', 'sleep_phase_5_min', 'timezone',
       'bedtime_start_delta', 'bedtime_end_delta', 'midpoint_at_delta',
       'heart_rate_5_min', 'hrv_5_min', 'readiness_temperature_deviation',
       'readiness_temperature_trend_deviation'],
      dtype='object')


# Columns to drop
# sleep_df
drop_col_sleep = ['efficiency', 'period', 'segment_state',
            'sleep_midpoint', 'sleep_phase_5_min', 'movement_30_sec',
            'timezone', 'bedtime_end_delta', 'midpoint_at_delta',
            'heart_rate_5_min', 'hrv_5_min'] # timezone? -- might need during EDA
sleep_df.drop(drop_col_sleep, axis = 1, inplace = True)

# biometric_df
drop_col_bio = ['average_met_minutes', 'equivalent_walking_distance', 
                'high_activity_met_minutes', 'inactivity_alerts', 
                'low_activity_met_minutes', 'medium_activity_met_minutes', 
                'sedentary_met_minutes', 'target_calories', 
                'target_meters', 'class_5_min', 
                'met_1_min', 'ring_met_1_min', 
                'temperature_trend_deviation']
biometric_df.drop(drop_col_bio, axis = 1, inplace = True)


# Reformat date columns
date_col_sleep = ['bedtime_end', 'bedtime_start', 'day']
# All times reformatted to UTC
sleep_df[date_col_sleep] = sleep_df[date_col_sleep].apply(pd.to_datetime, utc = True, errors = 'coerce')
biometric_df['day'] = biometric_df['day'].apply(pd.to_datetime, utc = True, errors = 'coerce')
daily_sleep_score_df['day'] = daily_sleep_score_df['day'].apply(pd.to_datetime, utc = True, errors = 'coerce')


print('biometric_df:')
print(biometric_df.select_dtypes(exclude = ['int64', 'float64']).dtypes)
print('\nsleep_df:')
print(sleep_df.select_dtypes(exclude = ['int64', 'float64']).dtypes)

biometric_df:
day    datetime64[ns, UTC]
dtype: object

sleep_df:
bedtime_end      datetime64[ns, UTC]
bedtime_start    datetime64[ns, UTC]
day              datetime64[ns, UTC]
type                          object
dtype: object


sleep_df['type'].unique()

array(['long_sleep', nan, 'late_nap', 'rest', 'sleep'], dtype=object)


# Initialize a time frame we can consider a nap/rest period (UTC format).
nap_upper = pd.to_datetime('23:59:00').time()
nap_lower = pd.to_datetime('15:00:00').time()
# Condition on the start, end time, and total sleep duration to filter out false-positives, remove any long_sleep types.
nap_bool = ((sleep_df['bedtime_start'].dt.time >= nap_lower) & 
            (sleep_df['bedtime_end'].dt.time <= nap_upper) & 
            (sleep_df['total_sleep_duration'] > 600) &
            (sleep_df['type'] != 'long_sleep'))
print('Naps identified:\n')
display(sleep_df[nap_bool])
# Insert a nap_today column and binary yes/no for each day.
# Initialize column with zeroes.
sleep_df['nap_today'] = np.zeros_like(sleep_df.shape[0])
# Use boolean array to identify nap days and iterate through to change nap_today to 1.
# TODO: This is not the pandas way but is fine with this amount of data. Use a better method later.
# Something like, sleep_df['nap_today'].where(~nap_bool, 1, inplace = True) does not work in this situation..
# Need each entry for that day to = 1 not just the nap itself.
nap_days = sleep_df[nap_bool]['day']
for day in nap_days:
    sleep_df.loc[sleep_df['day'] == day, 'nap_today'] = 1

Naps identified:


# Shift day one day into the future to align with relevant sleep score day.
biometric_df['day'] = biometric_df['day'] + pd.DateOffset(days = 1)
# Remove last day in biometric dataset since it won't have a sleep score with day shifting.
bio_max_idx = biometric_df.loc[biometric_df['day'] == biometric_df['day'].max()].index
biometric_df.drop(index = bio_max_idx, axis = 0, inplace = True)

# Verify each dataframe starts and ends on the same day.
print('Earliest Dates in each DF:')
print(biometric_df['day'].min(), '<-- biometric_df')
print(sleep_df['day'].min(), '<-- sleep_df')
print(daily_sleep_score_df['day'].min(), '<-- daily_sleep_score_df')
print('\nLatest Dates in each DF:')
print(biometric_df['day'].max(), '<-- biometric_df')
print(sleep_df['day'].max(), '<-- sleep_df')
print(daily_sleep_score_df['day'].max(), '<-- daily_sleep_score_df')

Earliest Dates in each DF:
2023-02-04 00:00:00+00:00 <-- biometric_df
2023-02-04 00:00:00+00:00 <-- sleep_df
2023-02-04 00:00:00+00:00 <-- daily_sleep_score_df

Latest Dates in each DF:
2024-02-06 00:00:00+00:00 <-- biometric_df
2024-02-06 00:00:00+00:00 <-- sleep_df
2024-02-06 00:00:00+00:00 <-- daily_sleep_score_df


# Identify days where there wasn't a long period of sleep
sleep_df_long = sleep_df.loc[sleep_df['type'] == 'long_sleep']

# Use merge() to find the differences in dataframes indicating days that didn't record a long sleep.
# Requires a method where the two df indexes don't match.
# Credit to: https://stackoverflow.com/questions/48647534/find-difference-between-two-data-frames
no_sleep = pd.DataFrame(sleep_df['day']).merge(sleep_df_long['day'], indicator = True, how='left').loc[lambda x : x['_merge']!='both']
print('Entries in which no long periods of sleep were recorded that day:')
display(no_sleep)

# We can remove these days with no long_sleep from the main dataframes as they won't have sleep scores anyways.
    # Sleep DF
sleep_df.drop(no_sleep.index, axis = 0, inplace = True)

    # Biometric DF
    # We can use our previous work and just remove the days that the two DFs don't share.
no_sleep_bio_days =  pd.DataFrame(biometric_df['day']).merge(sleep_df['day'], indicator = True, how='left').loc[lambda x : x['_merge']!='both']['day']
no_sleep_bio_index = biometric_df[biometric_df['day'].isin(no_sleep_bio_days)].index
biometric_df.drop(no_sleep_bio_index, axis = 0, inplace = True)

    # Daily Sleep Score DF
no_sleep_score_index = pd.DataFrame(daily_sleep_score_df['day']).merge(biometric_df['day'], indicator = True, how='left').loc[lambda x : x['_merge']!='both']['day'].index
daily_sleep_score_df.drop(index = no_sleep_score_index, axis = 0, inplace = True)

Entries in which no long periods of sleep were recorded that day:


# Instead of having 0-100 potential label predictions when using a classification model,
# we should create a categorical 'binned' label column.
# Create bins. We can use .5 to easily decide left edge or right edge issues.
bin_number = 10
bin_interval_even = np.linspace(0.5, 100.5, bin_number + 1, dtype = float)
bin_labels = [f"{i} to {i + bin_number}" for i in bin_interval_even if i != bin_interval_even[-1]]
# np.digitize() can map the sleep scores to the correct bins
bin_idx = np.digitize(daily_sleep_score_df['score'], bin_interval_even, right = False)
# Create new column with the new bin labels
bin_map_even = [bin_labels[value - 1] for value in bin_idx]
daily_sleep_score_df['score_bin_even'] = bin_map_even
print(bin_labels)
display(daily_sleep_score_df['score_bin_even'].value_counts())

daily_sleep_score_df['score'].plot(kind='kde')
plt.xlim(0, 100)
plt.vlines(bin_interval_even, ymin = 0, ymax = 0.055, color = 'purple', alpha = 0.5)
for idx, label in enumerate(bin_labels):
    plt.text(x = ((bin_interval_even[idx] + bin_interval_even[idx + 1]) / 2) - 1.5,
             y = (plt.ylim()[0] + plt.ylim()[1]) / 2,
             s = label,
             rotation = 'vertical',
             alpha = 0.5)
plt.xlabel('Sleep Score')

['0.5 to 10.5', '10.5 to 20.5', '20.5 to 30.5', '30.5 to 40.5', '40.5 to 50.5', '50.5 to 60.5', '60.5 to 70.5', '70.5 to 80.5', '80.5 to 90.5', '90.5 to 100.5']

score_bin_even
70.5 to 80.5     165
80.5 to 90.5     101
60.5 to 70.5      68
50.5 to 60.5       9
40.5 to 50.5       5
90.5 to 100.5      4
30.5 to 40.5       1
Name: count, dtype: int64

Text(0.5, 0, 'Sleep Score')


#labels_quantile = ['Worst', 'Bad', 'Below Avg','Above Avg', 'Good', 'Best']
labels_quantile = ['Worst', 'Bad', 'Avg', 'Good', 'Best']
bin_map_quantile, bin_interval_quantile = pd.qcut(daily_sleep_score_df['score'], 5, labels = labels_quantile, retbins = True)
daily_sleep_score_df['score_bin_quantile'] = bin_map_quantile

print('Quantile Bins:', bin_interval_quantile)
print('\nDistribution of Bins:', bin_map_quantile.value_counts())
daily_sleep_score_df['score'].plot(kind = 'kde')
plt.xlim(0, 100)
plt.vlines(bin_interval_quantile, ymin = 0, ymax = 0.055, color = 'purple', alpha = 0.5)
for idx, label in enumerate(labels_quantile):
    plt.text(x = ((bin_interval_quantile[idx] + bin_interval_quantile[idx + 1]) / 2) - 1.5,
             y = (plt.ylim()[0] + plt.ylim()[1]) / 2,
             s = f'{label} - {bin_map_quantile.value_counts()[label]}',
             rotation = 'vertical',
             alpha = 0.5)
plt.xlabel('Sleep Score')

Quantile Bins: [36. 69. 75. 78. 82. 94.]

Distribution of Bins: score
Bad      84
Worst    73
Best     68
Good     66
Avg      62
Name: count, dtype: int64

Text(0.5, 0, 'Sleep Score')


labels_custom = ['Worst', 'Bad', 'Below Avg', 'Above Avg', 'Good', 'Best']
median_score = daily_sleep_score_df['score'].median()
std_score = daily_sleep_score_df['score'].std()
bin_interval_custom = [0,
            median_score - (1.5 * std_score),
            median_score - (0.5 * std_score),
            median_score,
            median_score + (0.5 * std_score),
            median_score + (1.5 * std_score),
            100]
bin_map_custom = pd.cut(daily_sleep_score_df['score'], bins = bin_interval_custom, labels = labels_custom)
daily_sleep_score_df['score_bin_custom'] = bin_map_custom

print('Custom Quantile Bins:', bin_interval_custom)
print('\nDistribution of Bins:', bin_map_custom.value_counts())
daily_sleep_score_df['score'].plot(kind='kde')
plt.xlim(0, 100)
plt.vlines(bin_interval_custom, ymin = 0, ymax = 0.055, color = 'purple', alpha = 0.5)
for idx, label in enumerate(labels_custom):
    plt.text(x = ((bin_interval_custom[idx] + bin_interval_custom[idx + 1]) / 2) - 1.5,
             y = (plt.ylim()[0] + plt.ylim()[1]) / 2,
             s = f'{label} - {bin_map_custom.value_counts()[label]}',
             rotation = 'vertical',
             alpha = 0.5)
plt.xlabel('Sleep Score')

Custom Quantile Bins: [0, 64.39549639460978, 72.79849879820325, 77.0, 81.20150120179675, 89.60450360539022, 100]

Distribution of Bins: score
Below Avg    88
Good         78
Bad          77
Above Avg    75
Worst        30
Best          5
Name: count, dtype: int64

Text(0.5, 0, 'Sleep Score')


display(daily_sleep_score_df.tail(5))


# Sum each day's sleep metrics (important for days with interrupted sleep or naps).
# Overwrites each day's entry.
sleep_df = sleep_df[['day', 
         'deep_sleep_duration', 
         'light_sleep_duration', 
         'rem_sleep_duration', 
         'restless_periods', 
         'awake_time', 
         'time_in_bed', 
         'total_sleep_duration']].groupby(['day']).sum().reset_index()

# Features extracted only from the long sleep to include in the final dateframe.
sleep_df_long = sleep_df_long[['day', 
                               'average_breath', 
                               'average_heart_rate', 
                               'average_hrv', 
                               'latency', 
                               'lowest_heart_rate', 
                               'bedtime_start_delta',
                               'nap_today']]
# Merge the summed features with the long_sleep-only features.
sleep_df = sleep_df.merge(sleep_df_long, on = 'day', how = 'outer')
sleep_df.head(2)


print('Verify all DF.shape matches:', sleep_df.shape[0] == biometric_df.shape[0] == daily_sleep_score_df.shape[0])
print(sleep_df.shape, '<-- sleep_df')
print(biometric_df.shape, '<-- biometric_df')
print(daily_sleep_score_df.shape, '<-- daily_sleep_score_df')

Verify all DF.shape matches: True
(353, 15) <-- sleep_df
(353, 13) <-- biometric_df
(353, 5) <-- daily_sleep_score_df


# Merge sleep and biometric DFs.
bio_sleep_df = sleep_df.merge(biometric_df, on = 'day', how = 'outer')
# Merge sleep score df and feature df.
bio_sleep_df = daily_sleep_score_df.merge(bio_sleep_df, on = 'day', how = 'outer')
bio_sleep_df.reset_index(drop = True, inplace = True)

display(bio_sleep_df)
print(bio_sleep_df.columns)

Index(['score', 'day', 'score_bin_even', 'score_bin_quantile',
       'score_bin_custom', 'deep_sleep_duration', 'light_sleep_duration',
       'rem_sleep_duration', 'restless_periods', 'awake_time', 'time_in_bed',
       'total_sleep_duration', 'average_breath', 'average_heart_rate',
       'average_hrv', 'latency', 'lowest_heart_rate', 'bedtime_start_delta',
       'nap_today', 'active_calories', 'high_activity_time',
       'low_activity_time', 'medium_activity_time', 'meters_to_target',
       'non_wear_time', 'resting_time', 'sedentary_time', 'steps',
       'total_calories', 'temperature_deviation', 'spo2_percentage'],
      dtype='object')


# Explore and deal with NaN values.
print('Columns with NaNs present:', bio_sleep_df.columns[bio_sleep_df.isna().any()].to_list())
print(f"{'':<5}Total number of NaNs: {bio_sleep_df.isna().sum(axis = 1).sum()}")

Columns with NaNs present: ['temperature_deviation', 'spo2_percentage']
     Total number of NaNs: 48


# A good option to deal with each feature individually.
# >>> values = {{"Feature 1": .median(), "Feature 2": 0, "Feature 3: .mean()}}
# >>> df.fillna(value=values)
# Only feature with NaN currently is SPO2 and temp. deviation, using median for now.
# TODO: Come back and tune each feature.
for column in bio_sleep_df:
    if bio_sleep_df[column].dtype == 'float64' or bio_sleep_df[column].dtype == 'int64':
        bio_sleep_df[column].fillna(bio_sleep_df[column].median(), inplace = True)
    else:
        continue
print(f'Number of NaN values remaining after adjustments: {bio_sleep_df.isna().sum(axis = 1).sum()}')

Number of NaN values remaining after adjustments: 0


bio_sleep_df.describe()


pair_columns_sleep = ['score', 'score_bin_custom', 'deep_sleep_duration', 
                'light_sleep_duration', 'rem_sleep_duration', 
                'latency', 'bedtime_start_delta']
sns.pairplot(bio_sleep_df[pair_columns_sleep], hue = 'score_bin_custom')

<seaborn.axisgrid.PairGrid at 0x177bebc10>


pair_columns_sleep2 = ['score', 'score_bin_custom', 'restless_periods', 
                'awake_time', 'average_breath', 'average_hrv',
                 'nap_today']
sns.pairplot(bio_sleep_df[pair_columns_sleep2], hue = 'score_bin_custom')

<seaborn.axisgrid.PairGrid at 0x2896682d0>


pair_columns_active = ['score', 'score_bin_custom',
                'high_activity_time','low_activity_time', 'medium_activity_time',
                'resting_time', 'sedentary_time']
sns.pairplot(bio_sleep_df[pair_columns_active], hue = 'score_bin_custom')

<seaborn.axisgrid.PairGrid at 0x28d9dbed0>


# Function to transform inputs from seconds to hours.
def sec_to_hours(seconds):
    return seconds / 60 / 60

def reg_scatter_plot(feat_list):
    """Prints to output (n x 2) subplot array of scatter plots where x = feature in inputted list vs y = the target variable "score". 
    Can handle any number of input models.

    Parameters:
        features (list): Iterable of feature names from bio_sleep_df in the form of a string.
    
    Returns:
        None. Prints plot to output."""
    
    fig, ax = plt.subplots(int(len(feat_list)/2), 2, figsize = (6, 6), sharey = True)
    for i, feat in enumerate(feat_list):

        # Create regression line
        m, b = np.polyfit(bio_sleep_df[feat], bio_sleep_df['score'], 1)

        # i iterates through rows and columns with // and % operators.
        ax[i//2 , i%2].scatter(bio_sleep_df[feat], bio_sleep_df['score'], alpha = 0.5, s = 15)
        ax[i//2 , i%2].plot(bio_sleep_df[feat], m*bio_sleep_df[feat] + b, color = 'purple')
        ax[i//2 , i%2].set_title(feat)

    fig.suptitle('Sleep Score vs Types of Sleep')
    fig.supxlabel('Duration (seconds)')
    fig.supylabel('Sleep Score')
    fig.tight_layout()
    
    return None


print('Types of Sleep in Hours:')
sec_to_hours(bio_sleep_df[['deep_sleep_duration', 
                           'light_sleep_duration', 
                           'rem_sleep_duration', 
                           'total_sleep_duration']].describe())

Types of Sleep in Hours:


x_min = min(sec_to_hours(bio_sleep_df['total_sleep_duration']))
x_max = max(sec_to_hours(bio_sleep_df['total_sleep_duration']))
sec_to_hours(bio_sleep_df['total_sleep_duration']).plot(kind = 'kde')
plt.title('Distribution of Total Amount of Sleep per Day')
plt.xlabel('Hours of Sleep')
plt.xticks(np.arange(round(x_min - 1), round(x_max + 1), 0.5))
plt.xlim(round(x_min - 1), round(x_max + 1))
plt.grid(alpha = 0.4)


ax = sns.boxplot(x = bio_sleep_df['day'].dt.month, y = daily_sleep_score_df['score'])
ax.set(title = 'Sleep Score by Month Boxplot', xlabel = 'month')
plt.show()


ax = sns.boxplot(x = bio_sleep_df['day'].dt.month, y = bio_sleep_df['awake_time'])
ax.set(title = 'Awake Time by Month Boxplot', xlabel = 'month')
plt.show()


sleep_features = ['total_sleep_duration', 'deep_sleep_duration', 'rem_sleep_duration', 'light_sleep_duration']

reg_scatter_plot(sleep_features)


# Remove day column in preparation for model training.
bio_sleep_df.drop(['day'], axis = 1, inplace = True)


# Correlation
# Include score column but drop categorical bins to see correlation between response and predictors too.
corr_matrix = bio_sleep_df.drop(['score_bin_even', 'score_bin_quantile', 'score_bin_custom'], axis = 1).corr()
# Plot matrix.
fig, ax = plt.subplots(figsize = (16, 16))
sns.heatmap(corr_matrix, ax = ax, annot = True, fmt = '.1f')

<Axes: >


bio_sleep_df[['bedtime_start_delta']].describe()


plt.scatter(bio_sleep_df['bedtime_start_delta'], bio_sleep_df['score'], alpha = 0.5)
plt.xlabel('Bedtime Start Delta (seconds)')
plt.ylabel('Sleep Score')

Text(0, 0.5, 'Sleep Score')


late_bedtime_scores = bio_sleep_df[bio_sleep_df['bedtime_start_delta'] > 0]['score']
early_bedtime_scores = bio_sleep_df[bio_sleep_df['bedtime_start_delta'] < 0]['score']

# Sample mean checks.
print(f"The sample mean of sleep scores considered to have late bedtimes: {late_bedtime_scores.mean()}")
print(f"\nThe sample mean of sleep scores considered to have early bedtimes: {early_bedtime_scores.mean()}\n")

The sample mean of sleep scores considered to have late bedtimes: 72.30630630630631

The sample mean of sleep scores considered to have early bedtimes: 77.10743801652893


# Set level of significance.
alpha = 0.05


fig, ax = plt.subplots(2, 2, figsize = (6, 6))
sns.histplot(late_bedtime_scores, kde = True, ax = ax[0, 0])
sm.qqplot(late_bedtime_scores, fit = 'scale', line = '45', alpha = 0.5 , ax = ax[0, 1])
ax[0,0].set_title('Late Bedtime Score: Density')
ax[0,1].set_title('Late Bedtime Score: QQ Plot')
sns.histplot(early_bedtime_scores, kde = True, ax = ax[1, 0])
sm.qqplot(early_bedtime_scores, fit = 'scale', line = '45', alpha = 0.5 , ax = ax[1, 1])
ax[1,0].set_title('Early Bedtime Score: Density')
ax[1,1].set_title('Early Bedtime Score: QQ Plot')

fig.tight_layout()
plt.show()


# TODO: If data is updated the following is not setup if sample conditions change drastically.
# Let's test if our samples display normality with the shapiro-wilk test.
# Null hypothesis H_0: samples come from a normal distribution.
print('Early Bedtime:')
t_stat_sw, p_value_sw = stats.shapiro(early_bedtime_scores)
if p_value_sw > alpha:
    print(f"Shapiro-Wilk test (Normality): \n{'':<5}with a p-value of {p_value_sw:.4e}, \n{'':<5}we fail to reject the null hypothesis: samples are from a normally distributed population.")
    print("\nSamples DO come from a normal distribution.")
else:
    print(f"Shapiro-Wilk test (Normality): \n{'':<5}with a p-value of {p_value_sw:.4e}, \n{'':<5}we successfully reject the null hypothesis: samples are from a normally distributed population.")
    print("\nSamples do NOT come from a normal distribution.")

print('\nLate Bedtime:')
t_stat_sw, p_value_sw = stats.shapiro(late_bedtime_scores)
if p_value_sw > alpha:
    print(f"Shapiro-Wilk test (Normality): \n{'':<5}with a p-value of {p_value_sw:.4e}, \n{'':<5}we fail to reject the null hypothesis: samples are from a normally distributed population.")
    print("\nSamples DO come from a normal distribution.")
else:
    print(f"Shapiro-Wilk test (Normality): \n{'':<5}with a p-value of {p_value_sw:.4e}, \n{'':<5}we successfully reject the null hypothesis: samples are from a normally distributed population.")
    print("\nSamples do NOT come from a normal distribution.")

Early Bedtime:
Shapiro-Wilk test (Normality): 
     with a p-value of 1.0500e-08, 
     we successfully reject the null hypothesis: samples are from a normally distributed population.

Samples do NOT come from a normal distribution.

Late Bedtime:
Shapiro-Wilk test (Normality): 
     with a p-value of 2.2177e-06, 
     we successfully reject the null hypothesis: samples are from a normally distributed population.

Samples do NOT come from a normal distribution.


# Transformations
# Square root first
late_bedtime_scores = np.sqrt(max(late_bedtime_scores + 1) - late_bedtime_scores)
early_bedtime_scores = np.sqrt(max(early_bedtime_scores + 1) - early_bedtime_scores)


fig, ax = plt.subplots(2, 2, figsize = (6, 6))
sns.histplot(late_bedtime_scores, kde = True, ax = ax[0, 0])
sm.qqplot(late_bedtime_scores, fit = 'scale', line = '45', alpha = 0.5 , ax = ax[0, 1])
ax[0,0].set_title('Late Bedtime Score: Density')
ax[0,1].set_title('Late Bedtime Score: QQ Plot')
sns.histplot(early_bedtime_scores, kde = True, ax = ax[1, 0])
sm.qqplot(early_bedtime_scores, fit = 'scale', line = '45', alpha = 0.5 , ax = ax[1, 1])
ax[1,0].set_title('Early Bedtime Score: Density')
ax[1,1].set_title('Early Bedtime Score: QQ Plot')

fig.tight_layout()
plt.show()

# Null hypothesis H_0: samples come from a normal distribution.
print('Early Bedtime:')
t_stat_sw, p_value_sw = stats.shapiro(early_bedtime_scores)
if p_value_sw > alpha:
    print(f"Shapiro-Wilk test (Normality): \n{'':<5}with a p-value of {p_value_sw}, \n{'':<5}we fail to reject the null hypothesis: samples are from a normally distributed population.")
    print("\nSamples DO come from a normal distribution.")
else:
    print(f"Shapiro-Wilk test (Normality): \n{'':<5}with a p-value of {p_value_sw}, \n{'':<5}we successfully reject the null hypothesis: samples are from a normally distributed population.")
    print("\nSamples do NOT come from a normal distribution.")

print('\nLate Bedtime:')
t_stat_sw, p_value_sw = stats.shapiro(late_bedtime_scores)
if p_value_sw > alpha:
    print(f"Shapiro-Wilk test (Normality): \n{'':<5}with a p-value of {p_value_sw}, \n{'':<5}we fail to reject the null hypothesis: samples are from a normally distributed population.")
    print("\nSamples DO come from a normal distribution.")
else:
    print(f"Shapiro-Wilk test (Normality): \n{'':<5}with a p-value of {p_value_sw}, \n{'':<5}we successfully reject the null hypothesis: samples are from a normally distributed population.")
    print("\nSamples do NOT come from a normal distribution.")

Early Bedtime:
Shapiro-Wilk test (Normality): 
     with a p-value of 0.051339421421289444, 
     we fail to reject the null hypothesis: samples are from a normally distributed population.

Samples DO come from a normal distribution.

Late Bedtime:
Shapiro-Wilk test (Normality): 
     with a p-value of 0.06582079827785492, 
     we fail to reject the null hypothesis: samples are from a normally distributed population.

Samples DO come from a normal distribution.


# Calculate variance of each group.
s1 = np.var(late_bedtime_scores, ddof = 1) # ddof = 1 to calculate sample var.
s2 = np.var(early_bedtime_scores, ddof = 1)

# Null hypothesis H_0: Var(late) == Var(early)
# Larger var needs to be on top.
if s2 > s1:
    f_stat = s2 / s1
    df_num = len(early_bedtime_scores) - 1 # numerator
    df_denom = len(late_bedtime_scores) - 1 # denominator
else:
    f_stat = s1 / s2
    df_num = len(late_bedtime_scores) - 1 # numerator
    df_denom = len(early_bedtime_scores) - 1 # denominator

p_value_f = 1 - stats.f.cdf(f_stat, df_num, df_denom)
if p_value_f > alpha:
    print(f"F-test (variance): \n{'':<5}with a p-value of {p_value_f}, \n{'':<5}we fail to reject the null hypothesis: Var(late) == Var(early)")
    print("\nPopulations HAVE equal variances.")
else:
    print(f"F-test (variance): \n{'':<5}with a p-value of {p_value_f}, \n{'':<5}we successfully reject the null hypothesis: Var(late) == Var(early)")
    print("\nPopulations DO NOT have equal variances.")

F-test (variance): 
     with a p-value of 0.4487558615463687, 
     we fail to reject the null hypothesis: Var(late) == Var(early)

Populations HAVE equal variances.


# Finally let's perform a t-test to see if the mean scores are equal or not depending on the sleep time.
# Null hypothesis H_0: mean(late) == mean(early)
t_stat_t, p_value_t = stats.ttest_ind(late_bedtime_scores, early_bedtime_scores, alternative = 'two-sided')
if p_value_t > alpha:
    print(f"Two-Sample t-test (mean): \n{'':<5}with a p-value of {p_value_t}, \n{'':<5}we fail to reject the null hypothesis: mean(late) == mean(early)")
    print(f"\nThe means of the two samples are equal!")
else:
    print(f"Two-Sample t-test (mean): \n{'':<5}with a p-value of {p_value_t}, \n{'':<5}we successfully reject the null hypothesis: mean(late) == mean(early)")
    print(f"\nThe two sleep score population means are NOT equal!")

Two-Sample t-test (mean): 
     with a p-value of 0.013629025835902385, 
     we successfully reject the null hypothesis: mean(late) == mean(early)

The two sleep score population means are NOT equal!


linear_clf = sm.OLS(bio_sleep_df.iloc[:, :1], sm.add_constant(bio_sleep_df.iloc[:, 4:])).fit()
linear_clf.summary()


plt.scatter(linear_clf.fittedvalues, linear_clf.resid, alpha = 0.5)
plt.xlabel('Fitted Values')
plt.ylabel('Residuals')
plt.title('MLR Residuals vs Fitted Values')
plt.hlines(y = 0, xmin = 40, xmax = 100, colors = 'black', linestyles='--')
plt.show()


# Remove sign on correlation score for comparison.
# TODO: This method arbitrarily selects the first feature (of the two with a high score)
#   that it comes across. Change this method to have an insight-based selection method.
#   Possibly a preference list.
corr_matrix_abs = corr_matrix.abs()
# Drop score from x and y axis.
corr_matrix_abs.drop('score', axis = 0, inplace = True)
corr_matrix_abs.drop('score', axis = 1, inplace = True)
# Mask off only one of the correlation triangles to avoid duplicates.
corr_triangle = corr_matrix_abs.where(np.triu(np.ones(corr_matrix_abs.shape), k = 1).astype(bool))
# Create list of features with correlation higher than a threshold.
# Threshold set to > 0.75
to_drop = [column for column in corr_triangle.columns if any(corr_triangle[column] > 0.7)]
print(to_drop)

['total_sleep_duration', 'average_hrv', 'lowest_heart_rate', 'medium_activity_time', 'meters_to_target', 'steps', 'total_calories']


# For now let's manually remove features (see TODO note above).
choices_to_drop = [('average_hrv', 'average_heart_rate'),
           ('lowest_heart_rate', 'average_hrv'),
           ('medium_activity_time', 'active_calories'),
           ('meters_to_target', 'active_calories'),
           ('steps', 'a lot'),
           ('total_calories', 'a lot')]

bio_sleep_df.drop(['lowest_heart_rate', 
                   'active_calories', 
                   'meters_to_target', 
                   'steps', 
                   'total_calories', 
                   'average_heart_rate', 
                   'time_in_bed', 
                   'total_sleep_duration'],
                  axis = 1, inplace = True)


# Correlation
# Include score column to see correlation between response and predictors too.
corr_matrix = bio_sleep_df.drop(['score_bin_even', 'score_bin_quantile', 'score_bin_custom'], axis = 1).corr()
# Plot matrix.
fig, ax = plt.subplots(figsize = (16, 16))
sns.heatmap(corr_matrix, ax = ax, annot = True, fmt = '.1f')

<Axes: >


# TODO: Implement automated feature selection.


# Identify any score_bin_even where that class n == 1
display(bio_sleep_df['score_bin_even'].value_counts())
one_class = (bio_sleep_df['score_bin_even'].value_counts() == 1)
one_class_drop = (one_class.loc[one_class == True].index)

score_bin_even
70.5 to 80.5     165
80.5 to 90.5     101
60.5 to 70.5      68
50.5 to 60.5       9
40.5 to 50.5       5
90.5 to 100.5      4
30.5 to 40.5       1
Name: count, dtype: int64


# Remove outlier for now so we can class-balance the splits (min per class required = 2)
for value in one_class_drop:
    class_outlier_loc = bio_sleep_df.loc[bio_sleep_df['score_bin_even'] == value]
    bio_sleep_df.drop(class_outlier_loc.index, axis = 0, inplace = True)
    
bio_sleep_df.reset_index(drop = True, inplace = True)


rand_state = 87654321 #123456789 #1111 # Used in all models that accept a random state.
test_ratio = 0.2 
bin_type = 'score_bin_custom'

# Discretized Scores
# Stratify-balance the classes with the bin column (note this will slightly alter regression scores too).
X_train, X_test, y_train, y_test = train_test_split(bio_sleep_df.iloc[:, 4:], bio_sleep_df.iloc[:, :4], 
                                                    test_size = test_ratio, 
                                                    random_state = rand_state)#, 
                                                    #stratify = bio_sleep_df['score_bin_custom'])

# Binned/Categorical Class Labels
y_train_bin, y_test_bin = y_train[bin_type], y_test[bin_type]


# Centered and Standardized
scaler = StandardScaler()
X_train_scaled = pd.DataFrame(scaler.fit_transform(X_train), columns = scaler.get_feature_names_out())
X_test_scaled = pd.DataFrame(scaler.transform(X_test), columns = scaler.get_feature_names_out())
# Squeeze y's to put into a series, same format as train_test_split
y_train_scaled = pd.DataFrame(scaler.fit_transform(np.array(y_train['score']).reshape(-1, 1)), columns = ['score_scaled']).squeeze(axis = 1)
y_test_scaled = pd.DataFrame(scaler.transform(np.array(y_test['score']).reshape(-1, 1)), columns = ['score_scaled']).squeeze(axis = 1)


# Calculate the random chance baseline for our classification models.
# 1 / number of classes
random_base_class = 1 / len(bin_map_custom.unique())
print(f"The random chance baseline for classification models is: {random_base_class}")
# Calculate the Zero Rate baseline for our classification models.
zero_rate_base_class = y_train_bin.value_counts().max() / len(y_train_bin)
print(f"The Zero Rate accuracy baseline for classification models that always chooses the most populous class is: {zero_rate_base_class}")

The random chance baseline for classification models is: 0.16666666666666666
The Zero Rate accuracy baseline for classification models that always chooses the most populous class is: 0.2491103202846975


# Calculate the Zero Rate baseline for our classification models.
y_pred_baseline = np.full_like(y_train['score'], fill_value = y_train['score'].median(), dtype = float)
mae_baseline_regress = mean_absolute_error(y_train['score'], y_pred_baseline)
print(f"The baseline Mean Absolute Error (MAE) for regression models is: {mae_baseline_regress}")

The baseline Mean Absolute Error (MAE) for regression models is: 6.249110320284697


class MidpointNormalize(Normalize):
    def __init__(self, vmin=None, vmax=None, midpoint=None, clip=False):
        self.midpoint = midpoint
        Normalize.__init__(self, vmin, vmax, clip)

    def __call__(self, value, clip=None):
        x, y = [self.vmin, self.midpoint, self.vmax], [0, 0.5, 1]
        return np.ma.masked_array(np.interp(value, x, y))
    
def plotSearchGrid(grid):
    scores = [x for x in grid.cv_results_["mean_test_score"]]
    scores = np.array(scores).reshape(len(grid.param_grid['max_depth']), len(grid.param_grid['max_leaf_nodes']))

    plt.figure(figsize=(10, 8))
    plt.subplots_adjust(left=.2, right=0.95, bottom=0.15, top=0.95)
    plt.imshow(scores, interpolation='nearest', cmap=plt.cm.plasma,
               norm=MidpointNormalize(vmin=0.2, midpoint=0.92))
    plt.xlabel('max_leaf_nodes')
    plt.ylabel('max_depth')
    plt.colorbar()
    plt.xticks(np.arange(len(grid.param_grid['max_leaf_nodes'])), grid.param_grid['max_leaf_nodes'], rotation=90)
    plt.yticks(np.arange(len(grid.param_grid['max_depth'])), grid.param_grid['max_depth'])
    plt.title('Validation Score')
    plt.show()
    
def plotSearchGridSVM(grid):
    scores = [x for x in grid.cv_results_["mean_test_score"]]
    scores = np.array(scores).reshape(len(grid.param_grid['C']), len(grid.param_grid['gamma']))

    plt.figure(figsize=(10, 8))
    plt.subplots_adjust(left=.2, right=0.95, bottom=0.15, top=0.95)
    plt.imshow(scores, interpolation='nearest', cmap=plt.cm.plasma,
               norm=MidpointNormalize(vmin=0.2, midpoint=0.92))
    plt.xlabel('Gamma')
    plt.ylabel('C')
    plt.colorbar()
    plt.xticks(np.arange(len(grid.param_grid['gamma'])), grid.param_grid['gamma'], rotation=90)
    plt.yticks(np.arange(len(grid.param_grid['C'])), grid.param_grid['C'])
    plt.title('Validation Score')
    plt.show()


# Initialize list of fitted model objects.
models_regression = []
models_classification = []


mod_linear_OLS = sm.OLS(y_train['score'].reset_index(drop = True), sm.add_constant(X_train_scaled)).fit()
models_regression.append(mod_linear_OLS)
mod_linear_OLS.summary()


mod_log_regression = LogisticRegression(solver = 'lbfgs', multi_class = 'multinomial', max_iter = 30000, random_state = rand_state).fit(X_train, y_train_bin)
models_classification.append(mod_log_regression)

mod_lin_regression = LinearRegression().fit(X_train, y_train['score'])
models_regression.append(mod_lin_regression)

mod_ridge_regression = Ridge(max_iter = 10000, random_state = rand_state).fit(X_train, y_train['score'])
models_regression.append(mod_ridge_regression)

mod_ridge_classifier = RidgeClassifier(max_iter = 10000, random_state = rand_state).fit(X_train, y_train_bin)
models_classification.append(mod_ridge_classifier)

mod_lasso_regression = Lasso(max_iter = 10000, random_state = rand_state).fit(X_train, y_train['score'])
models_regression.append(mod_lasso_regression)

/opt/homebrew/lib/python3.11/site-packages/sklearn/linear_model/_logistic.py:460: ConvergenceWarning: lbfgs failed to converge (status=1):
STOP: TOTAL NO. of f AND g EVALUATIONS EXCEEDS LIMIT.

Increase the number of iterations (max_iter) or scale the data as shown in:
    https://scikit-learn.org/stable/modules/preprocessing.html
Please also refer to the documentation for alternative solver options:
    https://scikit-learn.org/stable/modules/linear_model.html#logistic-regression
  n_iter_i = _check_optimize_result(


params = {'C': np.logspace(6, 10, base = 2), 'gamma': np.logspace(-10, 0, base = 2)}
grid_svm_regression = GridSearchCV(SVR(kernel = 'rbf'), param_grid = params, cv = 3).fit(X_train_scaled, y_train['score'])
mod_svm_regression = grid_svm_regression.best_estimator_
models_regression.append(mod_svm_regression)

grid_svm_classifier = GridSearchCV(SVC(kernel = 'rbf', random_state = rand_state, probability = True), param_grid = params, cv = 3).fit(X_train_scaled, y_train_bin)
mod_svm_classifier = grid_svm_classifier.best_estimator_
models_classification.append(mod_svm_classifier)


plotSearchGridSVM(grid_svm_regression)

scores = [x for x in grid_svm_regression.cv_results_["mean_test_score"]]

max_index = scores.index(max(scores))
max_score = max(scores)
opt_C = grid_svm_regression.cv_results_['param_C'][max_index]
opt_Gamma = grid_svm_regression.cv_results_['param_gamma'][max_index]

print(f'GridSearch results indicate the optimized parameters are, \nC: {opt_C}, Gamma: {opt_Gamma}, Validation R2 Score: {max_score}')

# OR use the attribute .best_params_
grid_svm_regression.best_params_

GridSearch results indicate the optimized parameters are, 
C: 1024.0, Gamma: 0.001124953927697749, Validation R2 Score: 0.8797566276176538

{'C': 1024.0, 'gamma': 0.001124953927697749}


dt_classifier = DecisionTreeClassifier(max_depth = None, max_leaf_nodes = None, random_state = rand_state).fit(X_train, y_train_bin)
models_classification.append(dt_classifier)
print(dt_classifier.score(X_train, y_train_bin))
print(dt_classifier.score(X_test, y_test_bin))

1.0
0.30985915492957744


adaboost_classifier = AdaBoostClassifier(DecisionTreeClassifier(max_depth = 1),
                   n_estimators = 100,
                   random_state = rand_state,
                   learning_rate = 0.5
                   ).fit(X_train, y_train_bin)
models_classification.append(adaboost_classifier)
print(adaboost_classifier.score(X_test, y_test_bin))
results_adaboost = list(zip(adaboost_classifier.predict(X_test), y_test['score']))

0.23943661971830985


params = {'max_depth': np.linspace(1, 25, 25, dtype = int), 'max_leaf_nodes': np.linspace(2, 98, 25, dtype = int)}
grid_random_forest_regression = GridSearchCV(RandomForestRegressor(n_estimators = 250, random_state=rand_state), 
                    param_grid = params, 
                    cv = 3, 
                    #verbose = 3, 
                    n_jobs = -1
                    ).fit(X_train, y_train['score'])
mod_random_forest_regression = grid_random_forest_regression.best_estimator_
models_regression.append(mod_random_forest_regression)


params = {'max_depth': np.linspace(1, 25, 25, dtype = int), 'max_leaf_nodes': np.linspace(2, 98, 25, dtype = int)}
grid_random_forest_classifier = GridSearchCV(RandomForestClassifier(n_estimators = 250, random_state=rand_state), 
                    param_grid = params, 
                    cv = 3, 
                    #verbose = 3, 
                    n_jobs = -1
                    ).fit(X_train, y_train_bin)
mod_random_forest_classifier = grid_random_forest_classifier.best_estimator_
models_classification.append(mod_random_forest_classifier)


plotSearchGrid(grid_random_forest_regression)


feat_importance = grid_random_forest_regression.best_estimator_.feature_importances_
feat_imp_std = np.std([tree.feature_importances_ for tree in grid_random_forest_regression.best_estimator_.estimators_], axis = 0)
forest_importances = pd.Series(feat_importance, index=bio_sleep_df.columns[4:])

fig, ax = plt.subplots()
forest_importances.plot.bar(yerr = feat_imp_std, ax = ax)
ax.set_title("Feature Importance Regressor")
ax.set_ylabel("Mean Decrease in Impurity")
fig.tight_layout()


plotSearchGrid(grid_random_forest_classifier)


feat_importance = grid_random_forest_classifier.best_estimator_.feature_importances_
feat_imp_std = np.std([tree.feature_importances_ for tree in grid_random_forest_classifier.best_estimator_.estimators_], axis = 0)
forest_importances = pd.Series(feat_importance, index=bio_sleep_df.columns[4:])

fig, ax = plt.subplots()
forest_importances.plot.bar(yerr = feat_imp_std, ax = ax)
ax.set_title("Feature Importance - Classifier")
ax.set_ylabel("Mean Decrease in Impurity")
fig.tight_layout()


def scores_classification(models):
    """Calculates and returns accuracy score from fitted classifier models.

            Parameters:
                models (list): Iterable of fitted classifier model objects.
            
            Returns:
                Dictionary where model name is the key and it's accuracy score is the value. 
                    {'model name' : accuracy}"""
    # Create dictionary to return dataframe of scores.
    model_scores = {}
    print(f'Using target variable bin interval: {bin_type}')
    for i, mod in enumerate(models):
        # Check if model is the SVC(), uses the scaled features.
        if 'SVC' in str(mod):
            X_test_temp = X_test_scaled
        else:
            X_test_temp = X_test

        model_name = str(mod).partition('(')[0]
        accuracy_score = mod.score(X_test_temp, y_test_bin)
        f1score = f1_score(y_test_bin, mod.predict(X_test_temp), average = 'weighted')
        # Set scores to model in dictionary.
        model_scores[model_name] = {'Accuracy': accuracy_score,
                                    'F1_Score(class_weighted)': f1score}

    return model_scores

def ols_scores(model):
    """Calculates and returns regression metrics for statsmodels OLS() fitted models.

    Parameters:
        model (object): Fitted regressor model object.
    
    Returns:
        Nested dictionary where model name is the key and the model evaluation metrics are nested within. 
            {'model name' : {'R2': score,
                            'Adj_R2': score,
                            'RMSE': score,
                            'MAE': score}}"""

    y_pred = model.predict(sm.add_constant(X_test_scaled))
    r2 = model.rsquared
    adj_r2 = model.rsquared_adj
    rmse = mean_squared_error(y_test['score'], y_pred, squared = False)
    mae = mean_absolute_error(y_test['score'], y_pred)

    model_score = {'R2': r2,
                    'Adj_R2': adj_r2,
                    'RMSE': rmse,
                    'MAE': mae}

    return model_score

def scores_regression(models):
    """Calculates and returns R2, Adjusted R2, Root Mean Squared Error (RMSE), and Mean Absolute Error (MAE) from fitted regressor models.

        Parameters:
            models (list): Iterable of fitted regressor model objects.
        
        Returns:
            Nested dictionary where model name is the key and the model evaluation metrics are nested within. 
                {'model name' : {'R2': score,
                                'Adj_R2': score,
                                'RMSE': score,
                                'MAE': score}}"""
    # Create dictionary to return dataframe of scores.
    model_scores = {}

    for mod in models:
        model_name = str(mod).partition('(')[0]

        # Check if model is the SVR(), uses the scaled features.
        if 'SVR' in model_name:
            X_test_temp = X_test_scaled
        elif 'statsmodels' in model_name:
            model_scores['LinearRegression(scaled)'] = ols_scores(mod)
            continue
        else:
            X_test_temp = X_test

        y_pred = mod.predict(X_test_temp)
        r2 = mod.score(X_test_temp, y_test['score'])
        n , p = X_test_temp.shape # number of samples and predictors
        adj_r2 = 1 - (((1 - r2) * (n - 1)) / (n - p - 1))
        rmse = mean_squared_error(y_test['score'], y_pred, squared = False)
        mae = mean_absolute_error(y_test['score'], y_pred)
        # Set scores to model in dictionary.
        model_scores[model_name] = {'R2': r2,
                                    'Adj_R2': adj_r2,
                                    'RMSE': rmse,
                                    'MAE': mae}

    return model_scores

def plot_predictions_vs_score_regression(models):
    """Prints to output (n x 3) subplot array of each model's predicted sleep scores and observed (truth) sleep scores for each X_test, y_test pair, sorted by descending sleep score.
    (where n = number of model objects). Can handle any number of input models.

    Parameters:
        models (list): Iterable of fitted regressor model objects.
    
    Returns:
        None. Prints plot to output."""
    
    fig, ax = plt.subplots(int(np.ceil(len(models) / 3)), 3, figsize = (11, 7))
    for i, mod in enumerate(models):
        # Check if model is the SVR() or OLS(), uses the scaled features.
        if 'SVR' in str(mod):
            X_test_temp = X_test_scaled
        else:
            X_test_temp = X_test
            
        # Check if model is a statsmodel object, if so use their score conventions.
        if 'statsmodels' in str(mod):
            temp_scores = ols_scores(mod)
            y_pred = mod.predict(sm.add_constant(X_test_scaled))
            r2 = temp_scores['R2']
            adj_r2 = temp_scores['Adj_R2']
            rmse = temp_scores['RMSE']
            mae = temp_scores['MAE']
            model_name = 'LinearRegression(scaled)'
        else:
            y_pred = mod.predict(X_test_temp)
            r2 = mod.score(X_test_temp, y_test['score'])
            n , p = X_test_temp.shape # number of samples and predictors
            adj_r2 = 1 - (((1 - r2) * (n - 1)) / (n - p - 1))
            rmse = mean_squared_error(y_test['score'], y_pred, squared = False)
            mae = mean_absolute_error(y_test['score'], y_pred)
            model_name = str(mod).partition('(')[0]

        results = list(zip(y_pred, y_test['score']))
        results = sorted(results, key = lambda x: x[1], reverse = True)
        
        ax[i // 3, i % 3].plot(results, alpha = 0.85, drawstyle = 'steps')
        ax[i // 3, i % 3].set_xlabel('Index (sorted)')
        ax[i // 3, i % 3].set_ylabel('Sleep Score')
        ax[i // 3, i % 3].set_title(model_name)
        ax[i // 3, i % 3].grid(alpha = 0.4)
        ax[i // 3, i % 3].text(x = 0.1, 
                            y = 0.2, 
                            transform = ax[i // 3, i % 3].transAxes,
                            s = (f"{mae:0.3f} - MAE"
                            '\n'f"{rmse:0.3f} - RMSE"
                            '\n'f"{r2:0.3f} - R2"
                            '\n'f"{adj_r2:0.3f} - Adj_R2"))
        # Check if odd number of models and delete last subplot if true.
    if len(models) % 2 != 0:
        fig.delaxes(ax[len(models) // 3, (len(models) % 3)])
        
    plt.subplots_adjust(wspace = 0.2, hspace = 0.4)
    fig.legend(['Predicted', 'Truth'], loc = 'upper right')
    fig.suptitle('Predicted Score and True Score', fontsize = 12)

    return None

def plot_obs_vs_pred_regression(models):
    """Prints to output (n x 3) subplot array of each model's observed vs predicted plots.
    (where n = number of model objects). Can handle any number of input models.

    Parameters:
        models (list): Iterable of fitted regressor model objects.
    
    Returns:
        None. Prints plot to output."""
    
    fig, ax = plt.subplots(int(np.ceil(len(models) / 3)), 3, figsize = (11, 7))
    for i, mod in enumerate(models):
        # Check if model is the SVR(), uses the scaled features.
        if 'SVR' in str(mod):
            X_test_temp = X_test_scaled
        else:
            X_test_temp = X_test

        # Check if model is a statsmodel object, if so use their score conventions.
        if 'statsmodels' in str(mod):
            y_pred = mod.predict(sm.add_constant(X_test_scaled))
            model_name = 'LinearRegression(scaled)'
        else:
            y_pred = mod.predict(X_test_temp)
            model_name = str(mod).partition('(')[0]

        ax[i // 3, i % 3].scatter(y_pred, y_test['score'], alpha = 0.5, s = 15)
        ax[i // 3, i % 3].plot(np.arange(min(y_test['score']), 100, 1), np.arange(min(y_test['score']), 100, 1), 'k--', alpha = 0.6)
        ax[i // 3, i % 3].set_xlabel('Predictions')
        ax[i // 3, i % 3].set_ylabel('Observations')
        ax[i // 3, i % 3].set_title(model_name)
        ax[i // 3, i % 3].grid(alpha = 0.4)

        # Check if odd number of models and delete last subplot if true.
    if len(models) % 2 != 0:
        fig.delaxes(ax[len(models) // 3, (len(models) % 3)])
        
    plt.subplots_adjust(wspace = 0.3, hspace = 0.4)
    fig.suptitle('Observed (Truth) vs Predicted', fontsize = 12)

    return None

def plot_models_roc(models):
    """Prints to output (n x 3) subplot array of each model's ROC curve plot.
    (where n = number of model objects). Can handle any number of input models.

    Parameters:
        models (list): Iterable of fitted classifier model objects.
    
    Returns:
        None. Prints plot to output."""
    # Need to use LabelBinarizer() since we have multiple classes.
    label_binarizer = LabelBinarizer().fit(y_train_bin)
    y_onehot_test = label_binarizer.transform(y_test_bin)
    for mod in models:
        if 'Ridge' in str(mod): # Can't calculate probabilities in Ridge.
            continue
        if 'SVC' in str(mod): # Use scaled X_test values
            y_scores = mod.predict_proba(X_test_scaled)
        else:
            y_scores = mod.predict_proba(X_test)
        fpr, tpr, _ = roc_curve(y_onehot_test.ravel(), y_scores.ravel())
        auc_score = auc(fpr, tpr)
        plt.plot(fpr, tpr, label = f"{str(mod).partition('(')[0]} - AUC = {auc_score:0.2f}")

    plt.plot(np.arange(0, 1.1, 0.1), np.arange(0, 1.1, 0.1), 'k--')
    plt.title('ROC - Micro-Average OvR')
    plt.xlabel('False Positive Rate')
    plt.ylabel('True Positive Rate')
    plt.grid(alpha = 0.4)
    plt.legend(loc = 'lower right')
    
    return None

# Function to iterate through classification models and plot confusion matrix together.
def plot_models_confusion_matrix(models):
    """Prints to output (n x 3) subplot array of each model's confusion matrix.
    (where n = number of model objects). Can handle any number of input models.

    Parameters:
        models (list): Iterable of fitted classifier model objects.
    
    Returns:
        None. Prints plot to output."""
    fig, ax = plt.subplots(int(np.ceil(len(models) / 3)), 3, figsize = (10, 8), sharey = True)
    for i, mod in enumerate(models):
        # Check if model is the SVC(), uses the scaled features.
        if 'SVC' in str(mod):
            X_test_temp = X_test_scaled 
        else:
            X_test_temp = X_test

        cm = confusion_matrix(y_test_bin, mod.predict(X_test_temp), labels = y_test_bin.sort_values().unique())
        disp = ConfusionMatrixDisplay(confusion_matrix = cm, display_labels = y_test_bin.sort_values().unique())
        disp.plot(ax = ax[i // 3, i % 3], xticks_rotation = 90)
        disp.ax_.set_title(str(mod).partition('(')[0])
        disp.im_.colorbar.remove()

    # Check if odd number of models and delete last subplot if true.
    if len(models) % 2 != 0:
        fig.delaxes(ax[i // 3, (len(models) % 3)])
    fig.tight_layout()
    plt.subplots_adjust(wspace = 0.2, hspace = 0)
    fig.colorbar(disp.im_, ax = ax)
    fig.suptitle('Confusion Matrices of Classifier Models')
    
    return None


results_regression = scores_regression(models_regression)
display(pd.DataFrame(results_regression).T.sort_values(by = 'MAE'))


print(f"The baseline Mean Absolute Error (MAE) for regression models is: {mae_baseline_regress}")

The baseline Mean Absolute Error (MAE) for regression models is: 6.249110320284697


plot_predictions_vs_score_regression(models_regression)


plot_obs_vs_pred_regression(models_regression)


results_classification = scores_classification(models_classification)
display(pd.DataFrame(results_classification).T.sort_values(by = 'Accuracy', ascending = False))

Using target variable bin interval: score_bin_custom


print(f"The Zero Rate baseline for classification models that always chooses the most populous class is: {zero_rate_base_class}")

The Zero Rate baseline for classification models that always chooses the most populous class is: 0.2491103202846975


plot_models_confusion_matrix(models_classification)


plot_models_roc(models_classification)


print('Best Regression Model (by Mean Absolute Error (MAE)):') 
display(pd.DataFrame(results_regression).T.sort_values(by = 'MAE').iloc[0:1, :])
print('-----------------------------------------------------')
print('Best Classification Model (by Accuracy Score):')
display(pd.DataFrame(results_classification).T.sort_values(by = 'Accuracy', ascending = False).iloc[0:1, :])

Best Regression Model (by Mean Absolute Error (MAE)):

-----------------------------------------------------
Best Classification Model (by Accuracy Score):


# Exported via command line using:

# # HTML
# jupyter nbconvert Super_Sleep_Notebook.ipynb --to html
# jupyter nbconvert Super_Sleep_Notebook.ipynb --to html --HTMLExporter.theme=dark

# # PDF
# jupyter nbconvert Super_Sleep_Notebook.ipynb --to webpdf
# jupyter nbconvert Super_Sleep_Notebook.ipynb --to webpdf --HTMLExporter.theme=dark

	average_breath	average_heart_rate	average_hrv	awake_time	bedtime_end	bedtime_start	day	deep_sleep_duration	efficiency	latency	...	readiness_contributors_body_temperature	readiness_contributors_hrv_balance	readiness_contributors_previous_day_activity	readiness_contributors_previous_night	readiness_contributors_recovery_index	readiness_contributors_resting_heart_rate	readiness_contributors_sleep_balance	readiness_score	readiness_temperature_deviation	readiness_temperature_trend_deviation
0	13.625	63.25	77.0	4440.0	2023-02-04T07:08:22.000-06:00	2023-02-03T22:40:22.000-06:00	2023-02-04	4650.0	85.0	990.0	...	90.0	NaN	96.0	74.0	97.0	94.0	NaN	89.0	-0.38	NaN
1	15.125	92.49	19.0	6210.0	2023-02-05T09:37:28.000-06:00	2023-02-04T23:54:28.000-06:00	2023-02-05	4590.0	82.0	180.0	...	100.0	NaN	82.0	79.0	100.0	59.0	NaN	78.0	-0.04	0.15
2	15.125	87.67	NaN	1920.0	2023-02-05T20:10:57.000-06:00	2023-02-05T19:36:57.000-06:00	2023-02-06	30.0	6.0	0.0	...	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
3	15.125	82.50	37.0	720.0	2023-02-05T20:36:02.000-06:00	2023-02-05T20:20:02.000-06:00	2023-02-06	0.0	25.0	0.0	...	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN

	average_breath	average_heart_rate	average_hrv	awake_time	bedtime_end	bedtime_start	day	deep_sleep_duration	latency	light_sleep_duration	lowest_heart_rate	rem_sleep_duration	restless_periods	time_in_bed	total_sleep_duration	type	bedtime_start_delta	readiness_temperature_deviation	readiness_temperature_trend_deviation
31	13.500	74.47	34.0	6270.0	2023-02-20 21:16:13+00:00	2023-02-20 19:05:13+00:00	2023-02-20 00:00:00+00:00	420.0	570.0	840.0	67.0	330.0	25.0	7860	1590.0	NaN	47113	NaN	NaN
34	13.875	70.14	58.0	3540.0	2023-02-21 23:27:03+00:00	2023-02-21 21:59:03+00:00	2023-02-21 00:00:00+00:00	120.0	3600.0	1140.0	60.0	480.0	23.0	5280	1740.0	NaN	57543	NaN	NaN
44	14.750	90.50	25.0	3690.0	2023-02-27 01:11:23+00:00	2023-02-26 23:54:23+00:00	2023-02-27 00:00:00+00:00	300.0	2760.0	540.0	90.0	90.0	24.0	4620	930.0	NaN	60863	NaN	NaN
53	14.750	NaN	NaN	1140.0	2023-03-01 23:27:00+00:00	2023-03-01 22:51:00+00:00	2023-03-01 00:00:00+00:00	570.0	750.0	420.0	NaN	30.0	11.0	2160	1020.0	NaN	60660	NaN	NaN
63	14.625	71.67	NaN	930.0	2023-03-06 22:03:32+00:00	2023-03-06 21:30:32+00:00	2023-03-06 00:00:00+00:00	630.0	390.0	420.0	70.0	0.0	31.0	1980	1050.0	NaN	55832	NaN	NaN
88	13.750	65.15	75.0	1530.0	2023-03-17 21:33:33+00:00	2023-03-17 20:17:33+00:00	2023-03-17 00:00:00+00:00	90.0	450.0	2940.0	62.0	0.0	88.0	4560	3030.0	NaN	55053	NaN	NaN
89	13.750	68.33	63.0	2070.0	2023-03-17 23:52:33+00:00	2023-03-17 22:48:33+00:00	2023-03-18 00:00:00+00:00	540.0	1410.0	780.0	63.0	450.0	12.0	3840	1770.0	late_nap	64113	NaN	NaN
116	13.625	70.40	70.0	3030.0	2023-03-30 20:27:01+00:00	2023-03-30 19:26:01+00:00	2023-03-30 00:00:00+00:00	180.0	690.0	450.0	68.0	0.0	10.0	3660	630.0	NaN	51961	NaN	NaN
152	14.125	61.20	107.0	750.0	2023-04-23 21:10:49+00:00	2023-04-23 20:29:49+00:00	2023-04-23 00:00:00+00:00	1260.0	330.0	450.0	56.0	0.0	15.0	2460	1710.0	sleep	55789	NaN	NaN
180	14.875	71.20	80.0	1080.0	2023-05-10 18:45:42+00:00	2023-05-10 17:44:42+00:00	2023-05-10 00:00:00+00:00	810.0	570.0	1770.0	67.0	0.0	29.0	3660	2580.0	sleep	45882	NaN	NaN
197	13.500	66.75	117.0	1080.0	2023-05-16 23:12:57+00:00	2023-05-16 22:33:57+00:00	2023-05-17 00:00:00+00:00	60.0	1020.0	1200.0	63.0	0.0	19.0	2340	1260.0	late_nap	63237	NaN	NaN
210	13.500	73.29	81.0	1680.0	2023-05-20 23:05:22+00:00	2023-05-20 22:17:22+00:00	2023-05-21 00:00:00+00:00	0.0	1770.0	1200.0	69.0	0.0	17.0	2880	1200.0	late_nap	62242	NaN	NaN
233	14.500	70.20	64.0	1440.0	2023-06-03 23:11:00+00:00	2023-06-03 22:14:00+00:00	2023-06-04 00:00:00+00:00	1590.0	1320.0	390.0	64.0	0.0	13.0	3420	1980.0	NaN	-20760	NaN	NaN
251	14.250	59.73	123.0	1650.0	2023-06-11 21:39:31+00:00	2023-06-11 20:30:31+00:00	2023-06-11 00:00:00+00:00	900.0	600.0	1590.0	53.0	0.0	30.0	4140	2490.0	NaN	59431	NaN	NaN
258	14.125	75.67	73.0	570.0	2023-06-17 21:25:32+00:00	2023-06-17 20:52:32+00:00	2023-06-17 00:00:00+00:00	150.0	150.0	1260.0	72.0	0.0	18.0	1980	1410.0	sleep	60752	NaN	NaN
261	14.625	56.83	123.0	3090.0	2023-06-18 20:52:32+00:00	2023-06-18 19:36:32+00:00	2023-06-18 00:00:00+00:00	150.0	1770.0	1320.0	52.0	0.0	15.0	4560	1470.0	NaN	56192	NaN	-0.13
264	15.375	NaN	NaN	1980.0	2023-06-19 18:20:32+00:00	2023-06-19 17:33:02+00:00	2023-06-19 00:00:00+00:00	120.0	1920.0	750.0	NaN	0.0	29.0	2850	870.0	NaN	48782	NaN	-0.14
298	16.875	NaN	NaN	1141.0	2023-07-09 20:44:02+00:00	2023-07-09 20:07:01+00:00	2023-07-09 00:00:00+00:00	0.0	1140.0	1080.0	NaN	0.0	31.0	2221	1080.0	NaN	58021	NaN	-0.03
322	13.375	NaN	NaN	1050.0	2023-07-22 19:31:02+00:00	2023-07-22 19:02:02+00:00	2023-07-22 00:00:00+00:00	0.0	0.0	690.0	NaN	0.0	23.0	1740	690.0	NaN	54122	NaN	0.11
323	13.375	73.73	69.0	2070.0	2023-07-22 21:10:02+00:00	2023-07-22 19:36:02+00:00	2023-07-22 00:00:00+00:00	1170.0	1140.0	1800.0	70.0	600.0	72.0	5640	3570.0	NaN	56162	NaN	0.11
325	14.250	56.38	103.0	2070.0	2023-07-23 21:47:02+00:00	2023-07-23 20:21:32+00:00	2023-07-23 00:00:00+00:00	1200.0	1080.0	1860.0	53.0	0.0	19.0	5130	3060.0	NaN	58892	NaN	0.26
347	14.500	64.60	96.0	1260.0	2023-08-03 21:04:32+00:00	2023-08-03 20:16:02+00:00	2023-08-03 00:00:00+00:00	450.0	990.0	1200.0	62.0	0.0	3.0	2910	1650.0	NaN	58562	NaN	-0.14
349	14.750	72.75	79.0	1260.0	2023-08-04 22:19:30+00:00	2023-08-04 21:28:30+00:00	2023-08-05 00:00:00+00:00	300.0	600.0	1500.0	66.0	0.0	32.0	3060	1800.0	late_nap	62910	NaN	-0.02
353	14.125	NaN	NaN	691.0	2023-08-06 20:23:34+00:00	2023-08-06 19:54:33+00:00	2023-08-06 00:00:00+00:00	0.0	390.0	1050.0	NaN	0.0	10.0	1741	1050.0	NaN	57273	NaN	0.10
377	16.250	NaN	NaN	571.0	2023-08-20 22:38:00+00:00	2023-08-20 22:10:29+00:00	2023-08-20 00:00:00+00:00	90.0	420.0	960.0	NaN	30.0	0.0	1651	1080.0	NaN	61829	NaN	0.05
408	15.750	NaN	NaN	1202.0	2023-09-04 23:57:04+00:00	2023-09-04 23:25:32+00:00	2023-09-05 00:00:00+00:00	60.0	0.0	630.0	NaN	0.0	23.0	1892	690.0	NaN	-20068	NaN	0.16
452	17.875	89.67	21.0	1020.0	2023-09-25 16:37:57+00:00	2023-09-25 15:56:27+00:00	2023-09-26 00:00:00+00:00	210.0	690.0	1260.0	88.0	0.0	11.0	2490	1470.0	NaN	64587	NaN	0.06
455	15.375	79.80	NaN	1830.0	2023-09-27 16:11:57+00:00	2023-09-27 15:21:27+00:00	2023-09-28 00:00:00+00:00	30.0	0.0	1170.0	79.0	0.0	33.0	3030	1200.0	NaN	62487	NaN	0.18
463	16.375	NaN	NaN	570.0	2023-09-30 21:22:29+00:00	2023-09-30 20:46:29+00:00	2023-10-01 00:00:00+00:00	210.0	300.0	1380.0	NaN	0.0	8.0	2160	1590.0	NaN	-8011	NaN	0.14
508	13.375	61.17	101.0	420.0	2023-10-22 23:16:01+00:00	2023-10-22 22:30:01+00:00	2023-10-23 00:00:00+00:00	270.0	720.0	2010.0	58.0	60.0	22.0	2760	2340.0	late_nap	63001	NaN	-0.28
522	13.875	62.00	92.0	2670.0	2023-10-27 21:00:04+00:00	2023-10-27 19:52:34+00:00	2023-10-27 00:00:00+00:00	60.0	2130.0	1230.0	59.0	90.0	16.0	4050	1380.0	NaN	53554	NaN	0.03
540	15.375	68.50	66.0	1380.0	2023-11-05 22:15:01+00:00	2023-11-05 21:33:01+00:00	2023-11-05 00:00:00+00:00	360.0	1380.0	780.0	66.0	0.0	31.0	2520	1140.0	sleep	55981	NaN	-0.03
543	14.125	NaN	NaN	841.0	2023-11-06 23:03:30+00:00	2023-11-06 22:33:59+00:00	2023-11-06 00:00:00+00:00	30.0	810.0	900.0	NaN	0.0	7.0	1771	930.0	NaN	59639	NaN	NaN
590	15.500	90.60	NaN	1260.0	2023-11-30 20:52:33+00:00	2023-11-30 20:20:03+00:00	2023-11-30 00:00:00+00:00	0.0	0.0	690.0	86.0	0.0	23.0	1950	690.0	NaN	51603	NaN	NaN
594	14.375	73.23	51.0	6150.0	2023-12-02 23:08:02+00:00	2023-12-02 20:21:02+00:00	2023-12-02 00:00:00+00:00	1230.0	4860.0	1800.0	66.0	840.0	0.0	10020	3870.0	NaN	51662	0.47	0.53
617	14.750	74.90	45.0	1830.0	2023-12-16 22:11:30+00:00	2023-12-16 21:04:30+00:00	2023-12-16 00:00:00+00:00	150.0	1260.0	2040.0	66.0	0.0	38.0	4020	2190.0	NaN	54270	NaN	0.10
654	14.500	57.82	123.0	1500.0	2024-01-02 23:22:30+00:00	2024-01-02 21:40:30+00:00	2024-01-02 00:00:00+00:00	720.0	990.0	2910.0	54.0	990.0	16.0	6120	4620.0	sleep	56430	NaN	0.27
722	14.625	63.80	82.0	2190.0	2024-02-03 21:12:02+00:00	2024-02-03 20:07:02+00:00	2024-02-03 00:00:00+00:00	60.0	1470.0	1650.0	56.0	0.0	27.0	3900	1710.0	sleep	50822	NaN	-0.54
725	14.000	NaN	NaN	3060.0	2024-02-04 22:54:30+00:00	2024-02-04 21:45:00+00:00	2024-02-04 00:00:00+00:00	0.0	2880.0	1110.0	NaN	0.0	0.0	4170	1110.0	NaN	56700	NaN	-0.63

Score	Day (sleep_df)	Bedtime Start	Bedtime End	/	/	Day (biometric_df)	Steps	Active Calories
...	...	...	...	/	/	...	...	...
81	01/20/2024	1/19/2024 22:30	1/20/2024 07:00	/	/	01/20/2024	13000	800
65	01/21/2024	1/20/2024 23:45	1/21/2024 05:20	/	/	01/21/2024	21000	1900
...	...	...	...	/	/	...	...	...

Score	Day (sleep_df)	Bedtime Start	Bedtime End	/	/	~~Day (biometric_df)/~~	Steps	Active Calories
...	...	...	...	/	/	... ⇩	... ⇩	...⇩
81	01/20/2024	1/19/2024 22:30	1/20/2024 07:00	/	/	...⇩	...⇩	...⇩
65	01/21/2024	1/20/2024 23:45	1/21/2024 05:20	/	/	~~01/20/2024~~ +1 ⇩	13000 ⇩	800 ⇩
...	...	...	...	/	/	~~01/21/2024~~ +1 ⇩	21000 ⇩	1900 ⇩

	day	_merge
103	2023-03-24 00:00:00+00:00	left_only
140	2023-04-15 00:00:00+00:00	left_only
161	2023-05-01 00:00:00+00:00	left_only
162	2023-05-01 00:00:00+00:00	left_only
163	2023-05-02 00:00:00+00:00	left_only
173	2023-05-08 00:00:00+00:00	left_only
223	2023-05-28 00:00:00+00:00	left_only
224	2023-05-28 00:00:00+00:00	left_only
326	2023-07-24 00:00:00+00:00	left_only
402	2023-08-31 00:00:00+00:00	left_only
441	2023-09-21 00:00:00+00:00	left_only
442	2023-09-21 00:00:00+00:00	left_only
536	2023-11-05 00:00:00+00:00	left_only
537	2023-11-05 00:00:00+00:00	left_only
538	2023-11-05 00:00:00+00:00	left_only
539	2023-11-05 00:00:00+00:00	left_only
540	2023-11-05 00:00:00+00:00	left_only
588	2023-11-29 00:00:00+00:00	left_only

	day	active_calories	average_met_minutes	equivalent_walking_distance	high_activity_met_minutes	high_activity_time	inactivity_alerts	low_activity_met_minutes	low_activity_time	medium_activity_met_minutes	...	temperature_trend_deviation	contributors_activity_balance	contributors_hrv_balance	contributors_previous_day_activity	contributors_previous_night	contributors_recovery_index	contributors_resting_heart_rate	contributors_sleep_balance	contributors_body_temperature	spo2_percentage
0	2023-02-03	140	1.18750	1969	0	0	1	86	8460	8	...	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
1	2023-02-04	650	1.59375	10728	16	120	0	277	23220	172	...	NaN	NaN	NaN	96.0	74.0	97.0	94.0	NaN	90.0	98.523

	score	day	score_bin_even	score_bin_quantile	score_bin_custom
350	68	2024-02-02 00:00:00+00:00	60.5 to 70.5	Worst	Bad
351	76	2024-02-03 00:00:00+00:00	70.5 to 80.5	Avg	Below Avg
352	75	2024-02-04 00:00:00+00:00	70.5 to 80.5	Bad	Below Avg
353	83	2024-02-05 00:00:00+00:00	80.5 to 90.5	Best	Good
354	83	2024-02-06 00:00:00+00:00	80.5 to 90.5	Best	Good

	day	deep_sleep_duration	light_sleep_duration	rem_sleep_duration	restless_periods	awake_time	time_in_bed	total_sleep_duration	average_breath	average_heart_rate	average_hrv	latency	lowest_heart_rate	bedtime_start_delta	nap_today
0	2023-02-04 00:00:00+00:00	4650.0	15570.0	5820.0	282.0	4440.0	30480	26040.0	13.625	63.25	77.0	990.0	56.0	-4778	0
1	2023-02-05 00:00:00+00:00	4590.0	14970.0	9210.0	240.0	6210.0	34980	28770.0	15.125	92.49	19.0	180.0	75.0	-332	0

	score	deep_sleep_duration	light_sleep_duration	rem_sleep_duration	restless_periods	awake_time	time_in_bed	total_sleep_duration	average_breath	average_heart_rate	...	low_activity_time	medium_activity_time	meters_to_target	non_wear_time	resting_time	sedentary_time	steps	total_calories	temperature_deviation	spo2_percentage
count	353.000000	353.000000	353.000000	353.000000	353.000000	353.000000	353.000000	353.000000	353.000000	353.000000	...	353.000000	353.000000	353.000000	353.000000	353.000000	353.000000	353.000000	353.000000	353.000000	353.000000
mean	75.597734	5265.042493	15382.776204	5475.042493	244.045326	6445.603399	32571.694051	26122.861190	13.884561	65.803088	...	17263.682720	3128.838527	-1493.484419	2969.745042	31003.342776	31844.362606	10819.963173	2940.813031	0.007450	97.842986
std	8.403002	1667.128745	2500.618719	1456.503659	54.461654	2899.808704	4112.171365	3212.363992	0.481592	7.836005	...	4813.519199	1855.161430	4482.521299	6154.118003	4777.300153	6612.383545	3447.932402	247.672000	0.216197	0.851728
min	36.000000	210.000000	9150.000000	1110.000000	88.000000	1590.000000	19740.000000	14040.000000	12.750000	49.400000	...	900.000000	120.000000	-38700.000000	0.000000	8280.000000	5040.000000	1601.000000	2229.000000	-0.710000	92.072000
25%	71.000000	4320.000000	13710.000000	4590.000000	211.000000	4555.000000	30221.000000	24060.000000	13.625000	60.240000	...	13980.000000	2160.000000	-3700.000000	0.000000	28800.000000	27360.000000	8522.000000	2802.000000	-0.110000	97.472000
50%	77.000000	5130.000000	15300.000000	5490.000000	239.000000	5940.000000	32460.000000	26280.000000	13.875000	64.170000	...	16860.000000	2820.000000	-800.000000	360.000000	31380.000000	32640.000000	10388.000000	2917.000000	0.020000	98.027000
75%	81.000000	6120.000000	17040.000000	6390.000000	280.000000	7560.000000	35345.000000	28080.000000	14.125000	69.730000	...	19920.000000	3720.000000	1100.000000	3480.000000	33900.000000	36660.000000	12708.000000	3060.000000	0.120000	98.419000
max	94.000000	10200.000000	21690.000000	9420.000000	492.000000	27456.000000	51036.000000	35520.000000	15.875000	96.680000	...	36240.000000	24660.000000	8900.000000	51840.000000	49260.000000	50160.000000	23382.000000	5036.000000	1.420000	99.200000

	deep_sleep_duration	light_sleep_duration	rem_sleep_duration	total_sleep_duration
count	0.098056	0.098056	0.098056	0.098056
mean	1.462512	4.272993	1.520845	7.256350
std	0.463091	0.694616	0.404584	0.892323
min	0.058333	2.541667	0.308333	3.900000
25%	1.200000	3.808333	1.275000	6.683333
50%	1.425000	4.250000	1.525000	7.300000
75%	1.700000	4.733333	1.775000	7.800000
max	2.833333	6.025000	2.616667	9.866667

Dep. Variable:	score	R-squared:	0.885
Model:	OLS	Adj. R-squared:	0.877
Method:	Least Squares	F-statistic:	101.0
Date:	Fri, 01 Mar 2024	Prob (F-statistic):	1.30e-137
Time:	13:15:23	Log-Likelihood:	-869.51
No. Observations:	353	AIC:	1791.
Df Residuals:	327	BIC:	1892.
Df Model:	25
Covariance Type:	nonrobust

	coef	std err	t	P>\|t\|	[0.025	0.975]
const	4.3622	28.532	0.153	0.879	-51.768	60.493
deep_sleep_duration	-0.0005	0.001	-0.550	0.583	-0.002	0.001
light_sleep_duration	-0.0014	0.001	-1.458	0.146	-0.003	0.000
rem_sleep_duration	-5.281e-05	0.001	-0.056	0.955	-0.002	0.002
restless_periods	-0.0221	0.004	-5.402	0.000	-0.030	-0.014
awake_time	-0.0053	0.004	-1.426	0.155	-0.013	0.002
time_in_bed	0.0047	0.004	1.253	0.211	-0.003	0.012
total_sleep_duration	-0.0019	0.003	-0.692	0.489	-0.007	0.004
average_breath	0.3119	0.478	0.653	0.514	-0.628	1.252
average_heart_rate	-0.0980	0.071	-1.389	0.166	-0.237	0.041
average_hrv	-0.0070	0.014	-0.487	0.626	-0.035	0.021
latency	-0.0021	0.000	-12.681	0.000	-0.002	-0.002
lowest_heart_rate	0.0248	0.060	0.413	0.680	-0.093	0.143
bedtime_start_delta	-0.0006	6.11e-05	-9.205	0.000	-0.001	-0.000
nap_today	-0.6323	0.557	-1.134	0.257	-1.729	0.464
active_calories	-0.0025	0.010	-0.257	0.797	-0.022	0.017
high_activity_time	0.0005	0.001	0.329	0.742	-0.002	0.003
low_activity_time	-0.0001	0.000	-0.642	0.521	-0.001	0.000
medium_activity_time	-0.0002	0.000	-0.490	0.624	-0.001	0.001
meters_to_target	-4.253e-05	0.000	-0.416	0.678	-0.000	0.000
non_wear_time	7.165e-05	0.000	0.565	0.572	-0.000	0.000
resting_time	0.0001	0.000	1.042	0.298	-0.000	0.000
sedentary_time	6.74e-05	0.000	0.490	0.624	-0.000	0.000
steps	2.323e-05	0.000	0.195	0.845	-0.000	0.000
total_calories	0.0065	0.009	0.719	0.473	-0.011	0.024
temperature_deviation	0.3037	0.828	0.367	0.714	-1.325	1.932
spo2_percentage	0.1547	0.210	0.736	0.462	-0.259	0.568

Omnibus:	18.471	Durbin-Watson:	1.919
Prob(Omnibus):	0.000	Jarque-Bera (JB):	23.825
Skew:	-0.438	Prob(JB):	6.71e-06
Kurtosis:	3.924	Cond. No.	1.00e+16

Omnibus:	17.170	Durbin-Watson:	2.142
Prob(Omnibus):	0.000	Jarque-Bera (JB):	19.924
Skew:	-0.525	Prob(JB):	4.72e-05
Kurtosis:	3.775	Cond. No.	11.4

	R2	Adj_R2	RMSE	MAE
SVR	0.866523	0.820319	2.526703	2.041957
Lasso	0.838919	0.783160	2.775705	2.212797
Ridge	0.837268	0.780937	2.789897	2.220879
LinearRegression	0.837241	0.780901	2.790127	2.221564
LinearRegression(scaled)	0.881864	0.873748	2.790127	2.221564
RandomForestRegressor	0.665764	0.550067	3.998320	3.021399

	coef	std err	t	P>\|t\|	[0.025	0.975]
const	75.8826	0.178	425.234	0.000	75.531	76.234
deep_sleep_duration	3.6552	0.232	15.764	0.000	3.199	4.112
light_sleep_duration	3.4404	0.252	13.631	0.000	2.943	3.937
rem_sleep_duration	4.0485	0.192	21.048	0.000	3.670	4.427
restless_periods	-1.2579	0.247	-5.083	0.000	-1.745	-0.771
awake_time	-1.9036	0.210	-9.079	0.000	-2.316	-1.491
average_breath	0.0709	0.213	0.334	0.739	-0.348	0.489
average_hrv	0.3172	0.211	1.500	0.135	-0.099	0.734
latency	-2.5018	0.214	-11.704	0.000	-2.923	-2.081
bedtime_start_delta	-1.6799	0.224	-7.493	0.000	-2.121	-1.238
nap_today	-0.0134	0.194	-0.069	0.945	-0.395	0.368
high_activity_time	0.2556	0.196	1.307	0.192	-0.130	0.641
low_activity_time	0.0520	0.583	0.089	0.929	-1.097	1.201
medium_activity_time	0.5161	0.321	1.610	0.109	-0.115	1.147
non_wear_time	0.3917	0.675	0.580	0.562	-0.938	1.721
resting_time	0.6315	0.526	1.200	0.231	-0.405	1.668
sedentary_time	0.6850	0.809	0.847	0.398	-0.908	2.278
temperature_deviation	0.0768	0.197	0.390	0.697	-0.311	0.465
spo2_percentage	0.1773	0.202	0.878	0.381	-0.221	0.575

	Accuracy	F1_Score(class_weighted)
SVC	0.633803	0.626915
LogisticRegression	0.605634	0.607407
RandomForestClassifier	0.563380	0.538434
RidgeClassifier	0.436620	0.376509
DecisionTreeClassifier	0.309859	0.304995
AdaBoostClassifier	0.239437	0.174365

Sleep Quality Prediction - Supervised Machine Learning¶

Table of Contents ¶

1. Data Source Information ¶

2. Setup ¶

2.1. Environment Information for Reproducibility: ¶

2.2. Importing the Data: ¶

3. Data Cleaning and Processing ¶

3.1. First Looks: ¶

3.2. Remove Oura's Score Aggregates: ¶

3.3. Removing Irrelevant Features: ¶

3.4. Data Munging: ¶

3.5. Feature Engineering and Processing: ¶

3.5.1. Did You Take a Nap Today?: ¶

3.5.2. Aligning the Dataframes - Why workout today when there's always tomorrow?: ¶

3.5.3. Removing Days Without Sleep - No rest for the predicted: ¶

3.5.4. Binning a Target Variable?!: ¶

3.5.4.1. Even-Width Bins: ¶

3.5.4.2. Quantile/Median-Centered Bins: ¶

3.5.4.3. Custom Quantile Bins: ¶

3.5.5. Aggregating the Sleep Metrics: ¶

3.5.6. Joining the Dataframes: ¶

3.5.7. Dealing with Missing Data - NaNs: ¶

4. Exploratory Data Analysis (EDA) ¶

4.1. High-Level Summary: ¶

Pair Plots:¶

4.2. EDA Helper Functions: ¶

4.3. Visualizing the Sleep Metrics: ¶

4.4. Remove Day Column: ¶

4.5. Correlation Matrix: ¶

4.6. Analyzing Bedtime Delta - The moving target of sleep: ¶

4.6.1. Hypothesis Testing: ¶

4.6.1.1 Two-sample t-test Hypothesis: ¶

4.6.1.2 Checking for Normality: ¶

4.6.1.3 Checking for Variance Equivalence: ¶

4.6.1.4 Performing a Two-sample t-test: ¶

4.7. OLS Multiple Linear Regression as a Diagnostic Tool - Pre-Feature Removal: ¶

5. Feature Selection/Removal ¶

6. Outlier Removal ¶

7. Train Test Split ¶

7.1. Standardize Variables - Created separately on the same split: ¶

8. Supervised Learning Model Selection ¶

8.1. Model Baseline Evaluation Metrics: ¶

8.1.1. Classification Baseline: ¶

8.1.2. Regression Baseline: ¶

8.2. Model Helper Functions: ¶

8.3. Model Initializations: ¶

8.3.1. Linear/Logistic Regression: ¶

8.3.2. Support Vector Machine (SVM): ¶

8.3.3. Decision Trees: ¶

8.3.3.1 AdaBoost: ¶

8.3.3.2 Random Forest: ¶

9. Results ¶

9.1. Results Helper Functions: ¶

9.2. Results - Evaluation Metrics: ¶

9.2.1. Regression Models: ¶

9.2.2. Classification Models: ¶

10. Conclusions ¶

The Short Answer:¶

The Long Answer and Tangent:¶